//! GridGraph2D8C - 2D grid graph with 8-connectivity
//!
//! This module implements a 2D grid graph where each node connects to
//! its 8 neighbors (left, right, up, down, and 4 diagonals).
//!
//! 8-connectivity is useful for image segmentation where diagonal connections
//! are important, such as detecting diagonal edges or boundaries.

use crate::internal::memory::NodeIndex2D;
use crate::types::{GridCutError, Label, Parent8, Segment};
use num_traits::{Bounded, Num};
use std::collections::VecDeque;

/// 2D grid graph with 8-connected neighbors
///
/// This structure represents a regular 2D grid where each node can connect
/// to up to 8 neighbors (4 orthogonal + 4 diagonal). It implements the
/// Boykov-Kolmogorov maximum flow algorithm optimized for grid topology.
///
/// # Type Parameters
///
/// * `TCap` - Terminal capacity type (edge capacities to source/sink)
/// * `NCap` - Neighbor capacity type (edge capacities between grid nodes)
/// * `Flow` - Flow accumulation type (must be able to hold sums of capacities)
///
/// # Example
///
/// ```ignore
/// use grid_cut_rs::GridGraph2D8C;
///
/// let mut graph = GridGraph2D8C::<i32, i32, i32>::new(100, 100)?;
/// let node = graph.node_id(50, 50)?;
/// graph.set_terminal_cap(node, 10, 5)?;
/// graph.compute_maxflow();
/// let segment = graph.get_segment(node);
/// ```
#[derive(Debug)]
pub struct GridGraph2D8C<TCap, NCap, Flow>
where
    TCap: Num + Copy + PartialOrd,
    NCap: Num + Copy + PartialOrd,
    Flow: Num + Copy + PartialOrd + From<TCap> + From<NCap>,
{
    // Grid dimensions
    /// Original grid width (user-specified)
    width: usize,
    /// Original grid height (user-specified)
    height: usize,
    /// Padded width (rounded up to multiple of 8 for cache alignment)
    padded_width: usize,
    /// Padded height (rounded up to multiple of 8 for cache alignment)
    padded_height: usize,
    /// Number of 8x8 blocks in x direction
    width_blocks: usize,

    // Node indexing helper
    /// Helper for block-based node indexing and neighbor navigation
    indexer: NodeIndex2D,

    // Node state vectors (indexed by node_id)
    /// Node labels: Free, Source tree, or Sink tree
    labels: Vec<Label>,
    /// Parent edge direction in search tree
    parents: Vec<Parent8>,
    /// Parent node ID in search tree
    parent_ids: Vec<usize>,
    /// Timestamps for path validation (avoids redundant tree traversals)
    timestamps: Vec<u32>,

    // Edge capacity vectors
    /// Residual capacities for 8 neighbor directions
    /// rc[0] = GreaterX (right), rc[1] = GreaterY (down),
    /// rc[2] = LessY (up), rc[3] = LessX (left),
    /// rc[4] = GreaterXGreaterY (diagonal down-right), rc[5] = GreaterXLessY (diagonal up-right),
    /// rc[6] = LessXGreaterY (diagonal down-left), rc[7] = LessXLessY (diagonal up-left)
    rc: [Vec<NCap>; 8],
    /// Terminal capacities (source/sink edges)
    /// Positive = source capacity, Negative = sink capacity
    rc_st: Vec<TCap>,

    // Algorithm state
    /// Active node queue for BFS growth phase
    queue: VecDeque<usize>,
    /// Orphan nodes awaiting adoption
    orphans: VecDeque<usize>,
    /// Secondary orphan queue for two-level processing
    orphans2: VecDeque<usize>,
    /// Free nodes to reactivate after adoption
    free_nodes: VecDeque<usize>,

    /// Current timestamp for path validation
    time: u32,
    /// Accumulated maximum flow value
    max_flow: Flow,
}

impl<TCap, NCap, Flow> GridGraph2D8C<TCap, NCap, Flow>
where
    TCap: Num + Copy + PartialOrd + Bounded + std::ops::Neg<Output = TCap> + From<NCap>,
    NCap: Num + Copy + PartialOrd + Bounded + From<TCap>,
    Flow: Num + Copy + PartialOrd + From<TCap> + From<NCap>,
{
    /// Creates a new 2D 8-connected grid graph
    ///
    /// # Arguments
    ///
    /// * `width` - Grid width (must be > 0)
    /// * `height` - Grid height (must be > 0)
    ///
    /// # Returns
    ///
    /// Returns `Ok(GridGraph2D8C)` on success, or `Err(GridCutError)` if:
    /// - Width or height is 0
    /// - Memory allocation fails
    ///
    /// # Example
    ///
    /// ```ignore
    /// let graph = GridGraph2D8C::<i32, i32, i32>::new(640, 480)?;
    /// ```
    pub fn new(width: usize, height: usize) -> Result<Self, GridCutError> {
        // Validate dimensions
        if width == 0 || height == 0 {
            return Err(GridCutError::InvalidDimensions { width, height });
        }

        // Calculate padded dimensions (add 2 for boundary, then round to multiple of 8)
        let padded_width = crate::internal::memory::next_higher_mul8(width + 2);
        let padded_height = crate::internal::memory::next_higher_mul8(height + 2);
        let width_blocks = padded_width >> 3;

        // Create node indexer
        let indexer = NodeIndex2D::new(padded_width, padded_height);

        // Calculate total number of nodes (including padding)
        let num_nodes = padded_width * padded_height;

        // Initialize node state vectors
        let labels = vec![Label::Free; num_nodes];
        let parents = vec![Parent8::None; num_nodes];
        let parent_ids = vec![0; num_nodes];
        let timestamps = vec![0; num_nodes];

        // Initialize capacity vectors for 8 directions
        let rc = [
            vec![NCap::zero(); num_nodes], // GreaterX (right)
            vec![NCap::zero(); num_nodes], // GreaterY (down)
            vec![NCap::zero(); num_nodes], // LessY (up)
            vec![NCap::zero(); num_nodes], // LessX (left)
            vec![NCap::zero(); num_nodes], // GreaterXGreaterY (diagonal down-right)
            vec![NCap::zero(); num_nodes], // GreaterXLessY (diagonal up-right)
            vec![NCap::zero(); num_nodes], // LessXGreaterY (diagonal down-left)
            vec![NCap::zero(); num_nodes], // LessXLessY (diagonal up-left)
        ];
        let rc_st = vec![TCap::zero(); num_nodes];

        // Initialize algorithm state
        let queue = VecDeque::with_capacity(num_nodes);
        let orphans = VecDeque::with_capacity(num_nodes);
        let orphans2 = VecDeque::with_capacity(num_nodes);
        let free_nodes = VecDeque::with_capacity(num_nodes);

        Ok(Self {
            width,
            height,
            padded_width,
            padded_height,
            width_blocks,
            indexer,
            labels,
            parents,
            parent_ids,
            timestamps,
            rc,
            rc_st,
            queue,
            orphans,
            orphans2,
            free_nodes,
            time: 0,
            max_flow: Flow::zero(),
        })
    }

    /// Returns the node ID for grid coordinates (x, y)
    ///
    /// Coordinates are 0-based and must be within the grid bounds.
    ///
    /// # Arguments
    ///
    /// * `x` - X coordinate (0 to width-1)
    /// * `y` - Y coordinate (0 to height-1)
    ///
    /// # Returns
    ///
    /// Returns `Ok(node_id)` on success, or `Err(GridCutError::InvalidCoordinates)`
    /// if coordinates are out of bounds.
    ///
    /// # Example
    ///
    /// ```ignore
    /// let node = graph.node_id(10, 20)?;
    /// ```
    pub fn node_id(&self, x: usize, y: usize) -> Result<usize, GridCutError> {
        if x >= self.width || y >= self.height {
            return Err(GridCutError::InvalidCoordinates { x, y, z: None });
        }

        // Add 1 to account for boundary padding
        Ok(self.indexer.node_id(x + 1, y + 1))
    }

    /// Sets terminal capacities for a node
    ///
    /// Sets the edge capacities from the node to the source and sink terminals.
    /// This also initializes the node for processing by adding it to the free nodes queue.
    ///
    /// # Arguments
    ///
    /// * `node_id` - Node identifier (obtained from `node_id()`)
    /// * `cap_source` - Capacity of edge to source terminal
    /// * `cap_sink` - Capacity of edge to sink terminal
    ///
    /// # Returns
    ///
    /// Returns `Ok(())` on success, or `Err(GridCutError::InvalidNodeId)` if
    /// the node ID is invalid.
    ///
    /// # Example
    ///
    /// ```ignore
    /// let node = graph.node_id(10, 20)?;
    /// graph.set_terminal_cap(node, 100, 50)?;
    /// ```
    pub fn set_terminal_cap(
        &mut self,
        node_id: usize,
        cap_source: TCap,
        cap_sink: TCap,
    ) -> Result<(), GridCutError> {
        if node_id >= self.labels.len() {
            return Err(GridCutError::InvalidNodeId {
                node_id,
                max: self.labels.len() - 1,
            });
        }

        // Set terminal capacity (source - sink)
        self.rc_st[node_id] = cap_source - cap_sink;

        // Initialize node based on terminal capacity
        if cap_source > cap_sink {
            // More capacity to source - add to source tree
            self.labels[node_id] = Label::Source;
            self.parents[node_id] = Parent8::Terminal;
            self.parent_ids[node_id] = node_id;
            self.free_nodes.push_back(node_id);
        } else if cap_sink > cap_source {
            // More capacity to sink - add to sink tree
            self.labels[node_id] = Label::Sink;
            self.parents[node_id] = Parent8::Terminal;
            self.parent_ids[node_id] = node_id;
            self.free_nodes.push_back(node_id);
        }
        // If equal, node remains Free

        Ok(())
    }

    /// Sets neighbor edge capacity
    ///
    /// Sets the capacity of the edge from the node to a neighbor at the specified offset.
    ///
    /// # Arguments
    ///
    /// * `node_id` - Node identifier (obtained from `node_id()`)
    /// * `offset_x` - X offset to neighbor (-1, 0, or +1)
    /// * `offset_y` - Y offset to neighbor (-1, 0, or +1)
    /// * `cap` - Edge capacity
    ///
    /// # Valid Offsets for 8-connectivity
    ///
    /// * `(-1, 0)` - Left neighbor
    /// * `(+1, 0)` - Right neighbor
    /// * `(0, -1)` - Up neighbor
    /// * `(0, +1)` - Down neighbor
    /// * `(+1, +1)` - Diagonal down-right
    /// * `(+1, -1)` - Diagonal up-right
    /// * `(-1, +1)` - Diagonal down-left
    /// * `(-1, -1)` - Diagonal up-left
    ///
    /// # Returns
    ///
    /// Returns `Ok(())` on success, or an error if:
    /// - Node ID is invalid
    /// - Offset is not a valid 8-connected neighbor
    ///
    /// # Example
    ///
    /// ```ignore
    /// let node = graph.node_id(10, 20)?;
    /// graph.set_neighbor_cap(node, 1, 0, 50)?;  // Right neighbor
    /// graph.set_neighbor_cap(node, 1, 1, 30)?;  // Diagonal down-right
    /// ```
    pub fn set_neighbor_cap(
        &mut self,
        node_id: usize,
        offset_x: i32,
        offset_y: i32,
        cap: NCap,
    ) -> Result<(), GridCutError> {
        if node_id >= self.labels.len() {
            return Err(GridCutError::InvalidNodeId {
                node_id,
                max: self.labels.len() - 1,
            });
        }

        // Map offset to direction index
        let direction = match (offset_x, offset_y) {
            (1, 0) => 0,   // GreaterX (right)
            (0, 1) => 1,   // GreaterY (down)
            (0, -1) => 2,  // LessY (up)
            (-1, 0) => 3,  // LessX (left)
            (1, 1) => 4,   // GreaterXGreaterY (diagonal down-right)
            (1, -1) => 5,  // GreaterXLessY (diagonal up-right)
            (-1, 1) => 6,  // LessXGreaterY (diagonal down-left)
            (-1, -1) => 7, // LessXLessY (diagonal up-left)
            _ => {
                return Err(GridCutError::InvalidOffset {
                    offset_x,
                    offset_y,
                    offset_z: None,
                });
            }
        };

        self.rc[direction][node_id] = cap;
        Ok(())
    }

    /// Bulk set all capacities from arrays
    ///
    /// Sets all terminal and neighbor capacities in a single operation.
    /// This is more efficient than calling `set_terminal_cap` and `set_neighbor_cap`
    /// individually for each node.
    ///
    /// # Arguments
    ///
    /// * `cap_source` - Array of source terminal capacities (length = width * height)
    /// * `cap_sink` - Array of sink terminal capacities (length = width * height)
    /// * `cap_le` - Array of left edge capacities (length = width * height)
    /// * `cap_ge` - Array of right edge capacities (length = width * height)
    /// * `cap_el` - Array of up edge capacities (length = width * height)
    /// * `cap_eg` - Array of down edge capacities (length = width * height)
    /// * `cap_ne` - Array of diagonal up-right edge capacities (length = width * height)
    /// * `cap_se` - Array of diagonal down-right edge capacities (length = width * height)
    /// * `cap_sw` - Array of diagonal down-left edge capacities (length = width * height)
    /// * `cap_nw` - Array of diagonal up-left edge capacities (length = width * height)
    ///
    /// All arrays must be in row-major order: `index = y * width + x`
    ///
    /// # Returns
    ///
    /// Returns `Ok(())` on success, or `Err(GridCutError::InvalidCapacityArrayLength)`
    /// if any array has incorrect length.
    ///
    /// # Example
    ///
    /// ```ignore
    /// let size = width * height;
    /// let cap_source = vec![100; size];
    /// let cap_sink = vec![50; size];
    /// let cap_le = vec![10; size];
    /// let cap_ge = vec![10; size];
    /// let cap_el = vec![10; size];
    /// let cap_eg = vec![10; size];
    /// let cap_ne = vec![5; size];
    /// let cap_se = vec![5; size];
    /// let cap_sw = vec![5; size];
    /// let cap_nw = vec![5; size];
    ///
    /// graph.set_caps(&cap_source, &cap_sink, &cap_le, &cap_ge, &cap_el, &cap_eg,
    ///                &cap_ne, &cap_se, &cap_sw, &cap_nw)?;
    /// ```
    pub fn set_caps(
        &mut self,
        cap_source: &[TCap],
        cap_sink: &[TCap],
        cap_le: &[NCap],
        cap_ge: &[NCap],
        cap_el: &[NCap],
        cap_eg: &[NCap],
        cap_ne: &[NCap],
        cap_se: &[NCap],
        cap_sw: &[NCap],
        cap_nw: &[NCap],
    ) -> Result<(), GridCutError> {
        let expected_len = self.width * self.height;

        // Validate all array lengths
        for (_name, len) in [
            ("cap_source", cap_source.len()),
            ("cap_sink", cap_sink.len()),
            ("cap_le", cap_le.len()),
            ("cap_ge", cap_ge.len()),
            ("cap_el", cap_el.len()),
            ("cap_eg", cap_eg.len()),
            ("cap_ne", cap_ne.len()),
            ("cap_se", cap_se.len()),
            ("cap_sw", cap_sw.len()),
            ("cap_nw", cap_nw.len()),
        ] {
            if len != expected_len {
                return Err(GridCutError::InvalidCapacityArrayLength {
                    expected: expected_len,
                    actual: len,
                });
            }
        }

        // Set capacities for all nodes
        for y in 0..self.height {
            for x in 0..self.width {
                let idx = y * self.width + x;
                let node = self.indexer.node_id(x + 1, y + 1);

                // Set terminal capacity
                self.rc_st[node] = cap_source[idx] - cap_sink[idx];

                // Initialize node based on terminal capacity
                if cap_source[idx] > cap_sink[idx] {
                    self.labels[node] = Label::Source;
                    self.parents[node] = Parent8::Terminal;
                    self.parent_ids[node] = node;
                    self.free_nodes.push_back(node);
                } else if cap_sink[idx] > cap_source[idx] {
                    self.labels[node] = Label::Sink;
                    self.parents[node] = Parent8::Terminal;
                    self.parent_ids[node] = node;
                    self.free_nodes.push_back(node);
                }

                // Set neighbor capacities
                self.rc[0][node] = cap_ge[idx]; // GreaterX (right)
                self.rc[1][node] = cap_eg[idx]; // GreaterY (down)
                self.rc[2][node] = cap_el[idx]; // LessY (up)
                self.rc[3][node] = cap_le[idx]; // LessX (left)
                self.rc[4][node] = cap_se[idx]; // GreaterXGreaterY (diagonal down-right)
                self.rc[5][node] = cap_ne[idx]; // GreaterXLessY (diagonal up-right)
                self.rc[6][node] = cap_sw[idx]; // LessXGreaterY (diagonal down-left)
                self.rc[7][node] = cap_nw[idx]; // LessXLessY (diagonal up-left)
            }
        }

        Ok(())
    }

    /// Returns the segment (Source or Sink) that a node belongs to
    ///
    /// This method should be called after `compute_maxflow()` to determine
    /// which side of the cut each node is on.
    ///
    /// # Arguments
    ///
    /// * `node_id` - Node identifier (obtained from `node_id()`)
    ///
    /// # Returns
    ///
    /// Returns `Segment::Source` if the node is in the source segment,
    /// or `Segment::Sink` if the node is in the sink segment.
    ///
    /// # Panics
    ///
    /// Panics if node_id is out of bounds (this is a programming error).
    ///
    /// # Example
    ///
    /// ```ignore
    /// graph.compute_maxflow();
    /// let node = graph.node_id(10, 20)?;
    /// let segment = graph.get_segment(node);
    /// ```
    pub fn get_segment(&self, node_id: usize) -> Segment {
        match self.labels[node_id] {
            Label::Source => Segment::Source,
            Label::Sink => Segment::Sink,
            Label::Free => Segment::Sink, // Free nodes are assigned to sink by convention
        }
    }

    /// Returns the computed maximum flow value
    ///
    /// This method should be called after `compute_maxflow()` to get the
    /// total flow from source to sink.
    ///
    /// # Returns
    ///
    /// The maximum flow value.
    ///
    /// # Example
    ///
    /// ```ignore
    /// graph.compute_maxflow();
    /// let flow = graph.get_flow();
    /// println!("Maximum flow: {}", flow);
    /// ```
    pub fn get_flow(&self) -> Flow {
        self.max_flow
    }


    /// Computes the maximum flow from source to sink
    ///
    /// This method executes the Boykov-Kolmogorov maximum flow algorithm.
    /// After completion, use `get_flow()` to retrieve the flow value and
    /// `get_segment()` to determine which side of the cut each node is on.
    ///
    /// The algorithm consists of three phases that repeat until no augmenting
    /// path exists:
    /// 1. Growth: Search for an augmenting path using BFS
    /// 2. Augmentation: Push flow along the path
    /// 3. Adoption: Restore tree structure
    ///
    /// # Example
    ///
    /// ```ignore
    /// let mut graph = GridGraph2D8C::<i32, i32, i32>::new(100, 100)?;
    /// // ... set up capacities ...
    /// graph.compute_maxflow();
    /// println!("Max flow: {}", graph.get_flow());
    /// ```
    pub fn compute_maxflow(&mut self) {
        // Initialize active queue from free_nodes
        // Transfer all labeled nodes (Source or Sink) to the active queue
        while let Some(v) = self.free_nodes.pop_front() {
            self.queue.push_back(v);
        }

        // Main algorithm loop
        loop {
            // Growth phase: find augmenting path
            let path = match self.grow() {
                Some(p) => p,
                None => break, // No more augmenting paths - algorithm complete
            };

            // Increment timestamp for path validation
            self.time += 1;

            // Augmentation phase: push flow along path
            let (vs, vt, st) = path;
            self.augment(vs, vt, st);

            // Adoption phase: restore tree structure
            self.adopt();
        }

        // Final pass: classify any Free nodes by doing BFS from source
        // This ensures proper segmentation for the min-cut
        self.finalize_segmentation();
    }

    /// Final segmentation pass after maxflow completes
    /// The min-cut is determined by which nodes remain in the source tree
    /// Nodes that are Free are assigned to Sink by default
    fn finalize_segmentation(&mut self) {
        // Simply convert Free nodes to Sink
        // Nodes that are still labeled Source or Sink from the algorithm keep their labels
        for label in &mut self.labels {
            if *label == Label::Free {
                *label = Label::Sink;
            }
        }
    }

    // ========================================================================
    // Algorithm Implementation - Private Methods
    // ========================================================================

    /// Get neighbor node ID for a given direction
    #[inline]
    fn get_neighbor(&self, node: usize, direction: usize) -> usize {
        match direction {
            0 => self.indexer.neighbor_ge(node), // GreaterX (right)
            1 => self.indexer.neighbor_eg(node), // GreaterY (down)
            2 => self.indexer.neighbor_el(node), // LessY (up)
            3 => self.indexer.neighbor_le(node), // LessX (left)
            4 => {
                // GreaterXGreaterY (diagonal down-right): [+1, +1]
                let temp = self.indexer.neighbor_ge(node);
                self.indexer.neighbor_eg(temp)
            }
            5 => {
                // GreaterXLessY (diagonal up-right): [+1, -1]
                let temp = self.indexer.neighbor_ge(node);
                self.indexer.neighbor_el(temp)
            }
            6 => {
                // LessXGreaterY (diagonal down-left): [-1, +1]
                let temp = self.indexer.neighbor_le(node);
                self.indexer.neighbor_eg(temp)
            }
            7 => {
                // LessXLessY (diagonal up-left): [-1, -1]
                let temp = self.indexer.neighbor_le(node);
                self.indexer.neighbor_el(temp)
            }
            _ => unreachable!("Invalid direction"),
        }
    }

    /// Returns the sister (reverse) direction index
    #[inline]
    fn sister(direction: usize) -> usize {
        match direction {
            0 => 3, // GreaterX <-> LessX
            1 => 2, // GreaterY <-> LessY
            2 => 1, // LessY <-> GreaterY
            3 => 0, // LessX <-> GreaterX
            4 => 7, // GreaterXGreaterY <-> LessXLessY
            5 => 6, // GreaterXLessY <-> LessXGreaterY
            6 => 5, // LessXGreaterY <-> GreaterXLessY
            7 => 4, // LessXLessY <-> GreaterXGreaterY
            _ => unreachable!("Invalid direction"),
        }
    }

    /// Converts direction index to Parent8 enum
    #[inline]
    fn direction_to_parent(direction: usize) -> Parent8 {
        match direction {
            0 => Parent8::GreaterX,
            1 => Parent8::GreaterY,
            2 => Parent8::LessY,
            3 => Parent8::LessX,
            4 => Parent8::GreaterXGreaterY,
            5 => Parent8::GreaterXLessY,
            6 => Parent8::LessXGreaterY,
            7 => Parent8::LessXLessY,
            _ => unreachable!("Invalid direction"),
        }
    }

    /// Growth phase: searches for an augmenting path from source to sink
    ///
    /// Returns Some((source_node, sink_node, connecting_direction)) if a path is found,
    /// or None if no augmenting path exists.
    fn grow(&mut self) -> Option<(usize, usize, usize)> {
        while let Some(node) = self.queue.pop_front() {
            let label = self.labels[node];

            // Try all 8 neighbor directions
            for direction in 0..8 {
                let neighbor = self.get_neighbor(node, direction);

                if label == Label::Source {
                    // Expanding source tree
                    if self.rc[direction][node] > NCap::zero() {
                        match self.labels[neighbor] {
                            Label::Sink => {
                                // Found augmenting path!
                                return Some((node, neighbor, direction));
                            }
                            Label::Free => {
                                // Add free node to source tree
                                self.labels[neighbor] = Label::Source;
                                self.parents[neighbor] =
                                    Self::direction_to_parent(Self::sister(direction));
                                self.parent_ids[neighbor] = node;
                                self.queue.push_back(neighbor);
                            }
                            Label::Source => {
                                // Already in source tree, skip
                            }
                        }
                    }
                } else if label == Label::Sink {
                    // Expanding sink tree - check sister edge capacity
                    let sister_dir = Self::sister(direction);
                    if self.rc[sister_dir][neighbor] > NCap::zero() {
                        match self.labels[neighbor] {
                            Label::Source => {
                                // Found augmenting path!
                                return Some((neighbor, node, sister_dir));
                            }
                            Label::Free => {
                                // Add free node to sink tree
                                self.labels[neighbor] = Label::Sink;
                                self.parents[neighbor] =
                                    Self::direction_to_parent(Self::sister(direction));
                                self.parent_ids[neighbor] = node;
                                self.queue.push_back(neighbor);
                            }
                            Label::Sink => {
                                // Already in sink tree, skip
                            }
                        }
                    }
                }
            }
        }

        // No augmenting path found
        None
    }


    /// Find minimum residual capacity along source tree path
    fn find_minrf_s(&self, mut v: usize) -> NCap {
        let mut min_cap = NCap::max_value();

        while self.parents[v] != Parent8::Terminal {
            let parent = self.parents[v];
            let parent_id = self.parent_ids[v];

            // Get direction from parent to v
            let direction = match parent {
                Parent8::GreaterX => 0,
                Parent8::GreaterY => 1,
                Parent8::LessY => 2,
                Parent8::LessX => 3,
                Parent8::GreaterXGreaterY => 4,
                Parent8::GreaterXLessY => 5,
                Parent8::LessXGreaterY => 6,
                Parent8::LessXLessY => 7,
                _ => unreachable!("Invalid parent in source tree"),
            };

            let sister_dir = Self::sister(direction);
            let cap = self.rc[sister_dir][parent_id];

            if cap < min_cap {
                min_cap = cap;
            }

            v = parent_id;
        }

        // Check terminal capacity
        let terminal_cap = self.rc_st[v];
        if terminal_cap > TCap::zero() {
            let terminal_cap_n = NCap::from(terminal_cap);
            if terminal_cap_n < min_cap {
                min_cap = terminal_cap_n;
            }
        }

        min_cap
    }

    /// Find minimum residual capacity along sink tree path
    fn find_minrf_t(&self, mut v: usize) -> NCap {
        let mut min_cap = NCap::max_value();

        while self.parents[v] != Parent8::Terminal {
            let parent = self.parents[v];
            let parent_id = self.parent_ids[v];

            // Get direction from parent to v
            let direction = match parent {
                Parent8::GreaterX => 0,
                Parent8::GreaterY => 1,
                Parent8::LessY => 2,
                Parent8::LessX => 3,
                Parent8::GreaterXGreaterY => 4,
                Parent8::GreaterXLessY => 5,
                Parent8::LessXGreaterY => 6,
                Parent8::LessXLessY => 7,
                _ => unreachable!("Invalid parent in sink tree"),
            };

            let sister_dir = Self::sister(direction);
            let cap = self.rc[sister_dir][parent_id];

            if cap < min_cap {
                min_cap = cap;
            }

            v = parent_id;
        }

        // Check terminal capacity (negative for sink)
        let terminal_cap = self.rc_st[v];
        if terminal_cap < TCap::zero() {
            let terminal_cap_n = NCap::from(-terminal_cap);
            if terminal_cap_n < min_cap {
                min_cap = terminal_cap_n;
            }
        }

        min_cap
    }

    /// Augment flow along source tree path
    fn aug_s(&mut self, mut v: usize, bottleneck: NCap) {
        while self.parents[v] != Parent8::Terminal {
            let parent = self.parents[v];
            let parent_id = self.parent_ids[v];

            // Get direction from parent to v
            let direction = match parent {
                Parent8::GreaterX => 0,
                Parent8::GreaterY => 1,
                Parent8::LessY => 2,
                Parent8::LessX => 3,
                Parent8::GreaterXGreaterY => 4,
                Parent8::GreaterXLessY => 5,
                Parent8::LessXGreaterY => 6,
                Parent8::LessXLessY => 7,
                _ => unreachable!("Invalid parent in source tree"),
            };

            let sister_dir = Self::sister(direction);

            // Update residual capacities
            self.rc[sister_dir][parent_id] = self.rc[sister_dir][parent_id] - bottleneck;
            self.rc[direction][v] = self.rc[direction][v] + bottleneck;

            // If edge saturated, mark node as orphan
            if self.rc[sister_dir][parent_id] == NCap::zero() {
                self.parents[v] = Parent8::None;
                self.orphans.push_back(v);
            }

            v = parent_id;
        }

        // Update terminal capacity
        let bottleneck_t = TCap::from(bottleneck);
        self.rc_st[v] = self.rc_st[v] - bottleneck_t;

        // If terminal edge saturated, mark node as orphan
        if self.rc_st[v] == TCap::zero() {
            self.parents[v] = Parent8::None;
            self.orphans.push_back(v);
        }
    }

    /// Augment flow along sink tree path
    fn aug_t(&mut self, mut v: usize, bottleneck: NCap) {
        while self.parents[v] != Parent8::Terminal {
            let parent = self.parents[v];
            let parent_id = self.parent_ids[v];

            // Get direction from parent to v
            let direction = match parent {
                Parent8::GreaterX => 0,
                Parent8::GreaterY => 1,
                Parent8::LessY => 2,
                Parent8::LessX => 3,
                Parent8::GreaterXGreaterY => 4,
                Parent8::GreaterXLessY => 5,
                Parent8::LessXGreaterY => 6,
                Parent8::LessXLessY => 7,
                _ => unreachable!("Invalid parent in sink tree"),
            };

            let sister_dir = Self::sister(direction);

            // Update residual capacities
            self.rc[sister_dir][parent_id] = self.rc[sister_dir][parent_id] - bottleneck;
            self.rc[direction][v] = self.rc[direction][v] + bottleneck;

            // If edge saturated, mark node as orphan
            if self.rc[sister_dir][parent_id] == NCap::zero() {
                self.parents[v] = Parent8::None;
                self.orphans.push_back(v);
            }

            v = parent_id;
        }

        // Update terminal capacity (negative for sink)
        let bottleneck_t = TCap::from(bottleneck);
        self.rc_st[v] = self.rc_st[v] + bottleneck_t;

        // If terminal edge saturated, mark node as orphan
        if self.rc_st[v] == TCap::zero() {
            self.parents[v] = Parent8::None;
            self.orphans.push_back(v);
        }
    }

    /// Augmentation phase: push flow along the augmenting path
    fn augment(&mut self, vs: usize, vt: usize, st: usize) {
        // Find minimum residual capacity (bottleneck)
        let mut bottleneck = self.rc[st][vs];

        let min_s = self.find_minrf_s(vs);
        if min_s < bottleneck {
            bottleneck = min_s;
        }

        let min_t = self.find_minrf_t(vt);
        if min_t < bottleneck {
            bottleneck = min_t;
        }

        // Update edge connecting the two trees
        self.rc[st][vs] = self.rc[st][vs] - bottleneck;
        let sister_st = Self::sister(st);
        self.rc[sister_st][vt] = self.rc[sister_st][vt] + bottleneck;

        // Augment along source tree path
        self.aug_s(vs, bottleneck);

        // Augment along sink tree path
        self.aug_t(vt, bottleneck);

        // Accumulate flow
        self.max_flow = self.max_flow + Flow::from(bottleneck);
    }

    /// Find origin (terminal) of a tree path using timestamps
    /// Returns Parent8::Terminal if path leads to terminal, Parent8::None otherwise
    fn find_origin(&mut self, start: usize) -> Parent8 {
        let mut v = start;

        // Follow parent pointers until we reach a terminal or detect invalid path
        let mut visited = Vec::new();
        while self.parents[v] != Parent8::Terminal && self.parents[v] != Parent8::None {
            // Check if this node was already validated in current timestamp
            if self.timestamps[v] == self.time {
                // Mark all nodes on this path
                for &node in &visited {
                    self.timestamps[node] = self.time;
                }
                return Parent8::Terminal;
            }

            visited.push(v);
            v = self.parent_ids[v];

            // Safety check to prevent infinite loops
            if visited.len() > self.labels.len() {
                return Parent8::None;
            }
        }

        // Mark all nodes on this path with current timestamp if valid
        if self.parents[v] == Parent8::Terminal {
            for &node in &visited {
                self.timestamps[node] = self.time;
            }
            self.timestamps[v] = self.time;
        }

        self.parents[v]
    }

    /// Adoption phase: restore tree structure after augmentation
    fn adopt(&mut self) {
        self.orphans2.clear();
        self.free_nodes.clear();

        while !self.orphans.is_empty() || !self.orphans2.is_empty() {
            let v = if !self.orphans2.is_empty() {
                self.orphans2.pop_front().unwrap()
            } else {
                self.orphans.pop_back().unwrap()
            };

            let label = self.labels[v];

            // Try to find a new valid parent
            let mut found_parent = false;

            // First, check if this node can reconnect directly to terminal
            let has_terminal_capacity = if label == Label::Source {
                self.rc_st[v] > TCap::zero()
            } else {
                self.rc_st[v] < TCap::zero()
            };

            if has_terminal_capacity {
                // Can reconnect to terminal
                self.parents[v] = Parent8::Terminal;
                self.parent_ids[v] = v;
                self.timestamps[v] = self.time;
                found_parent = true;
            } else {
                // Try to find a neighbor parent
                for direction in 0..8 {
                    let neighbor = self.get_neighbor(v, direction);

                    // Check if neighbor can be a valid parent
                    if self.labels[neighbor] == label {
                        let sister_dir = Self::sister(direction);

                        // Check if edge from neighbor to v has capacity
                        if self.rc[sister_dir][neighbor] > NCap::zero() {
                            // Check if neighbor has valid path to terminal
                            if self.find_origin(neighbor) == Parent8::Terminal {
                                // Found valid parent!
                                self.parents[v] = Self::direction_to_parent(direction);
                                self.parent_ids[v] = neighbor;
                                self.timestamps[v] = self.time;
                                found_parent = true;
                                break;
                            }
                        }
                    }
                }
            }

            if !found_parent {
                // Could not find valid parent - make node free
                self.labels[v] = Label::Free;
                self.free_nodes.push_back(v);

                // Add neighbors that had v as parent to orphan queue
                for direction in 0..8 {
                    let neighbor = self.get_neighbor(v, direction);

                    if self.labels[neighbor] == label {
                        let sister_dir = Self::sister(direction);
                        let expected_parent = Self::direction_to_parent(sister_dir);

                        if self.parents[neighbor] == expected_parent
                            && self.parent_ids[neighbor] == v
                        {
                            self.parents[neighbor] = Parent8::None;
                            self.orphans2.push_back(neighbor);
                        }
                    }
                }
            }
        }

        // Re-activate neighbors of freed nodes
        // When a node becomes free, its neighbors in the trees might now be able to expand
        while let Some(v) = self.free_nodes.pop_front() {
            for direction in 0..8 {
                let neighbor = self.get_neighbor(v, direction);
                let label = self.labels[neighbor];

                // Only reactivate nodes that are in a tree (not free, not orphan)
                if (label == Label::Source || label == Label::Sink)
                    && self.parents[neighbor] != Parent8::None
                {
                    // Check if this neighbor can now expand into the freed node v
                    let can_expand = if label == Label::Source {
                        // Source tree: check if neighbor -> v edge has capacity
                        self.rc[direction][neighbor] > NCap::zero()
                    } else {
                        // Sink tree: check if v -> neighbor edge has capacity (sister direction)
                        let sister_dir = Self::sister(direction);
                        self.rc[sister_dir][neighbor] > NCap::zero()
                    };

                    if can_expand {
                        self.queue.push_back(neighbor);
                    }
                }
            }
        }
    }
}

// Implement GridGraph trait for GridGraph2D8C
impl<TCap, NCap, Flow> crate::types::GridGraph for GridGraph2D8C<TCap, NCap, Flow>
where
    TCap: Num + Copy + PartialOrd + Bounded + std::ops::Neg<Output = TCap> + From<NCap>,
    NCap: Num + Copy + PartialOrd + Bounded + From<TCap>,
    Flow: Num + Copy + PartialOrd + From<TCap> + From<NCap>,
{
    type TerminalCap = TCap;
    type NeighborCap = NCap;
    type Flow = Flow;

    fn compute_maxflow(&mut self) {
        self.compute_maxflow()
    }

    fn get_flow(&self) -> Self::Flow {
        self.get_flow()
    }

    fn get_segment(&self, node_id: usize) -> Segment {
        self.get_segment(node_id)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_constructor_valid_dimensions() {
        let graph = GridGraph2D8C::<i32, i32, i32>::new(10, 20);
        assert!(graph.is_ok());

        let graph = graph.unwrap();
        assert_eq!(graph.width, 10);
        assert_eq!(graph.height, 20);
        assert_eq!(graph.max_flow, 0);
    }

    #[test]
    fn test_constructor_zero_dimensions() {
        let result = GridGraph2D8C::<i32, i32, i32>::new(0, 10);
        assert!(result.is_err());

        let result = GridGraph2D8C::<i32, i32, i32>::new(10, 0);
        assert!(result.is_err());
    }

    #[test]
    fn test_node_id_valid_coordinates() {
        let graph = GridGraph2D8C::<i32, i32, i32>::new(10, 20).unwrap();

        assert!(graph.node_id(0, 0).is_ok());
        assert!(graph.node_id(9, 19).is_ok());
        assert!(graph.node_id(5, 10).is_ok());
    }

    #[test]
    fn test_node_id_invalid_coordinates() {
        let graph = GridGraph2D8C::<i32, i32, i32>::new(10, 20).unwrap();

        assert!(graph.node_id(10, 5).is_err());
        assert!(graph.node_id(5, 20).is_err());
    }

    #[test]
    fn test_set_terminal_cap() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(10, 10).unwrap();
        let node = graph.node_id(5, 5).unwrap();

        assert!(graph.set_terminal_cap(node, 100, 50).is_ok());
        assert_eq!(graph.labels[node], Label::Source);

        let node2 = graph.node_id(3, 3).unwrap();
        graph.set_terminal_cap(node2, 30, 80).unwrap();
        assert_eq!(graph.labels[node2], Label::Sink);
    }

    #[test]
    fn test_set_neighbor_cap_orthogonal() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(10, 10).unwrap();
        let node = graph.node_id(5, 5).unwrap();

        // Test orthogonal directions
        assert!(graph.set_neighbor_cap(node, 1, 0, 10).is_ok()); // Right
        assert!(graph.set_neighbor_cap(node, -1, 0, 20).is_ok()); // Left
        assert!(graph.set_neighbor_cap(node, 0, 1, 30).is_ok()); // Down
        assert!(graph.set_neighbor_cap(node, 0, -1, 40).is_ok()); // Up

        assert_eq!(graph.rc[0][node], 10); // GreaterX
        assert_eq!(graph.rc[3][node], 20); // LessX
        assert_eq!(graph.rc[1][node], 30); // GreaterY
        assert_eq!(graph.rc[2][node], 40); // LessY
    }

    #[test]
    fn test_set_neighbor_cap_diagonal() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(10, 10).unwrap();
        let node = graph.node_id(5, 5).unwrap();

        // Test diagonal directions
        assert!(graph.set_neighbor_cap(node, 1, 1, 15).is_ok()); // Down-right
        assert!(graph.set_neighbor_cap(node, 1, -1, 25).is_ok()); // Up-right
        assert!(graph.set_neighbor_cap(node, -1, 1, 35).is_ok()); // Down-left
        assert!(graph.set_neighbor_cap(node, -1, -1, 45).is_ok()); // Up-left

        assert_eq!(graph.rc[4][node], 15); // GreaterXGreaterY
        assert_eq!(graph.rc[5][node], 25); // GreaterXLessY
        assert_eq!(graph.rc[6][node], 35); // LessXGreaterY
        assert_eq!(graph.rc[7][node], 45); // LessXLessY
    }

    #[test]
    fn test_set_neighbor_cap_invalid_offset() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(10, 10).unwrap();
        let node = graph.node_id(5, 5).unwrap();

        // Invalid offset (too far)
        let result = graph.set_neighbor_cap(node, 2, 0, 10);
        assert!(result.is_err());

        let result = graph.set_neighbor_cap(node, 0, 2, 10);
        assert!(result.is_err());
    }

    #[test]
    fn test_set_caps_valid() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(3, 2).unwrap();

        let size = 3 * 2;
        let cap_source = vec![100; size];
        let cap_sink = vec![50; size];
        let cap_le = vec![10; size];
        let cap_ge = vec![20; size];
        let cap_el = vec![30; size];
        let cap_eg = vec![40; size];
        let cap_ne = vec![5; size];
        let cap_se = vec![6; size];
        let cap_sw = vec![7; size];
        let cap_nw = vec![8; size];

        let result = graph.set_caps(
            &cap_source, &cap_sink, &cap_le, &cap_ge, &cap_el, &cap_eg, &cap_ne, &cap_se,
            &cap_sw, &cap_nw,
        );
        assert!(result.is_ok());

        let node = graph.node_id(1, 1).unwrap();
        assert_eq!(graph.labels[node], Label::Source);
        assert_eq!(graph.rc[0][node], 20); // GreaterX
        assert_eq!(graph.rc[4][node], 6); // GreaterXGreaterY (SE)
    }

    #[test]
    fn test_set_caps_wrong_length() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(3, 2).unwrap();

        let wrong_size = 5;
        let cap_source = vec![100; wrong_size];
        let cap_sink = vec![50; 6];
        let cap_le = vec![10; 6];
        let cap_ge = vec![20; 6];
        let cap_el = vec![30; 6];
        let cap_eg = vec![40; 6];
        let cap_ne = vec![5; 6];
        let cap_se = vec![6; 6];
        let cap_sw = vec![7; 6];
        let cap_nw = vec![8; 6];

        let result = graph.set_caps(
            &cap_source, &cap_sink, &cap_le, &cap_ge, &cap_el, &cap_eg, &cap_ne, &cap_se,
            &cap_sw, &cap_nw,
        );
        assert!(result.is_err());
    }

    #[test]
    fn test_get_segment() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(10, 10).unwrap();
        let node = graph.node_id(5, 5).unwrap();

        graph.labels[node] = Label::Source;
        assert_eq!(graph.get_segment(node), Segment::Source);

        graph.labels[node] = Label::Sink;
        assert_eq!(graph.get_segment(node), Segment::Sink);

        graph.labels[node] = Label::Free;
        assert_eq!(graph.get_segment(node), Segment::Sink);
    }

    #[test]
    fn test_get_flow() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(10, 10).unwrap();
        assert_eq!(graph.get_flow(), 0);

        graph.max_flow = 42;
        assert_eq!(graph.get_flow(), 42);
    }

    // ========================================================================
    // Integration Tests for Maxflow Algorithm with 8-connectivity
    // ========================================================================

    #[test]
    fn test_maxflow_simple_2x2_with_diagonals() {
        // Test a 2x2 grid with diagonal connections
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(2, 2).unwrap();

        let n00 = graph.node_id(0, 0).unwrap();
        let n10 = graph.node_id(1, 0).unwrap();
        let n01 = graph.node_id(0, 1).unwrap();
        let n11 = graph.node_id(1, 1).unwrap();

        // Set terminal capacities
        graph.set_terminal_cap(n00, 100, 0).unwrap();
        graph.set_terminal_cap(n11, 0, 100).unwrap();

        // Set orthogonal capacities
        graph.set_neighbor_cap(n00, 1, 0, 5).unwrap(); // [0,0] -> [1,0]
        graph.set_neighbor_cap(n00, 0, 1, 5).unwrap(); // [0,0] -> [0,1]
        graph.set_neighbor_cap(n10, 0, 1, 5).unwrap(); // [1,0] -> [1,1]
        graph.set_neighbor_cap(n01, 1, 0, 5).unwrap(); // [0,1] -> [1,1]

        // Add diagonal connection
        graph.set_neighbor_cap(n00, 1, 1, 10).unwrap(); // [0,0] -> [1,1] diagonal

        // Compute maxflow
        graph.compute_maxflow();

        // Flow should be positive (multiple paths available)
        let flow = graph.get_flow();
        assert!(flow > 0, "Flow should be positive");

        // Source and sink should be in different segments
        assert_eq!(graph.get_segment(n00), Segment::Source);
        assert_eq!(graph.get_segment(n11), Segment::Sink);
    }

    #[test]
    fn test_maxflow_diagonal_only_path() {
        // Test where the only path is through diagonals
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(3, 3).unwrap();

        let n00 = graph.node_id(0, 0).unwrap();
        let n11 = graph.node_id(1, 1).unwrap();
        let n22 = graph.node_id(2, 2).unwrap();

        // Set terminal capacities
        graph.set_terminal_cap(n00, 30, 0).unwrap();
        graph.set_terminal_cap(n22, 0, 30).unwrap();

        // Only set diagonal connections (no orthogonal)
        graph.set_neighbor_cap(n00, 1, 1, 15).unwrap(); // [0,0] -> [1,1]
        graph.set_neighbor_cap(n11, 1, 1, 15).unwrap(); // [1,1] -> [2,2]

        // Compute maxflow
        graph.compute_maxflow();

        // Flow should be limited by diagonal capacities
        let flow = graph.get_flow();
        assert_eq!(flow, 15);

        assert_eq!(graph.get_segment(n00), Segment::Source);
        assert_eq!(graph.get_segment(n22), Segment::Sink);
    }

    #[test]
    fn test_maxflow_compare_4c_vs_8c() {
        // Compare 4-connected and 8-connected on same problem
        // 8-connected should find more flow due to diagonal paths

        // First, test with 4-connectivity
        let mut graph_4c = crate::grid_cut::grid_2d_4c::GridGraph2D4C::<i32, i32, i32>::new(3, 3).unwrap();

        let n00_4c = graph_4c.node_id(0, 0).unwrap();
        let n22_4c = graph_4c.node_id(2, 2).unwrap();

        graph_4c.set_terminal_cap(n00_4c, 100, 0).unwrap();
        graph_4c.set_terminal_cap(n22_4c, 0, 100).unwrap();

        // Set uniform orthogonal capacities
        for y in 0..3 {
            for x in 0..3 {
                let node = graph_4c.node_id(x, y).unwrap();
                if x < 2 {
                    graph_4c.set_neighbor_cap(node, 1, 0, 10).unwrap();
                }
                if y < 2 {
                    graph_4c.set_neighbor_cap(node, 0, 1, 10).unwrap();
                }
            }
        }

        graph_4c.compute_maxflow();
        let flow_4c = graph_4c.get_flow();

        // Now test with 8-connectivity
        let mut graph_8c = GridGraph2D8C::<i32, i32, i32>::new(3, 3).unwrap();

        let n00_8c = graph_8c.node_id(0, 0).unwrap();
        let n22_8c = graph_8c.node_id(2, 2).unwrap();

        graph_8c.set_terminal_cap(n00_8c, 100, 0).unwrap();
        graph_8c.set_terminal_cap(n22_8c, 0, 100).unwrap();

        // Set same orthogonal capacities plus diagonals
        for y in 0..3 {
            for x in 0..3 {
                let node = graph_8c.node_id(x, y).unwrap();
                if x < 2 {
                    graph_8c.set_neighbor_cap(node, 1, 0, 10).unwrap();
                }
                if y < 2 {
                    graph_8c.set_neighbor_cap(node, 0, 1, 10).unwrap();
                }
                // Add diagonal connections
                if x < 2 && y < 2 {
                    graph_8c.set_neighbor_cap(node, 1, 1, 10).unwrap();
                }
            }
        }

        graph_8c.compute_maxflow();
        let flow_8c = graph_8c.get_flow();

        // 8-connected should have at least as much flow as 4-connected
        assert!(flow_8c >= flow_4c);
        
        // Both should have positive flow
        assert!(flow_4c > 0);
        assert!(flow_8c > 0);
    }

    #[test]
    fn test_maxflow_no_flow_disconnected() {
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(3, 3).unwrap();

        let source = graph.node_id(0, 0).unwrap();
        let sink = graph.node_id(2, 2).unwrap();

        graph.set_terminal_cap(source, 100, 0).unwrap();
        graph.set_terminal_cap(sink, 0, 100).unwrap();

        // No neighbor capacities - disconnected

        graph.compute_maxflow();
        assert_eq!(graph.get_flow(), 0);
    }

    #[test]
    fn test_maxflow_capacity_constraints() {
        // Verify all residual capacities are non-negative after maxflow
        let mut graph = GridGraph2D8C::<i32, i32, i32>::new(2, 2).unwrap();

        let n00 = graph.node_id(0, 0).unwrap();
        let n11 = graph.node_id(1, 1).unwrap();

        graph.set_terminal_cap(n00, 15, 0).unwrap();
        graph.set_terminal_cap(n11, 0, 15).unwrap();

        // Set all capacities
        graph.set_neighbor_cap(n00, 1, 0, 10).unwrap();
        graph.set_neighbor_cap(n00, 0, 1, 10).unwrap();
        graph.set_neighbor_cap(n00, 1, 1, 10).unwrap();

        graph.compute_maxflow();

        // All residual capacities should be non-negative
        for direction in 0..8 {
            for node_cap in &graph.rc[direction] {
                assert!(*node_cap >= 0);
            }
        }
    }
}
